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ABSTRACT 


A comiprehensive computer code has been developed through the use of the 
governing equations of motion and the Boussinesq hypothesis and the rise of a vortex 
pair in a linearly density-stratified medium has been calculated through the use of the 
Biot-Savart law and a finite difference scheme. The results have shown that the code ts 
capable of predicting the behavior of the vortices and the propagation of the internal 
waves resulting from the vortex motion. The reasons for the finite mse of the vortices 
in a Stratified medium is clearly explained in terms of the counter-rotating vortices 


induced by the stratification. 
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I. INTRODUCTION 


Vortices and vortex wakes have become a major theme of aerodynamics research 


since the advent of the large aircraft and the understanding of their evolution required 


an eXamination of many fundamental problems in fluid mechanics. Much of the 


progress made during the past two decades was discussed at the Symposium on 


Aircraft Wake Turbulence and Its Detection [Ref. 1] and at the Aircraft Wake Vortices 


Conference [Ref. 2]. Comprehensive reviews of the entire subject have been given by 
Donaldson and Bilanin [Ref. 3], Widnall (Ref. 4] and Hallock and Eberle (Ref. 5]. 


These studies, as well as numerous others carried out since 1977, have uncovered 


a number of complex problems which must be resolved in order to achieve a better 


understanding of the important features of -trailing vortices in homogeneous and 


stratified media. The principal ones are as follows [Ref. 6]: 


(a) 


(b) 


(¢) 


(d) 


(e) 


Roll-up process: The velocity and turbulence distribution at anv station 
behind the wing depend onthe wing section. wing-tip shape, Revnolds 
number, wing incidence. and the distancé of the station from the wing (Ref. 7]. 

he distributions of the initial velocity and turbulence, which influence the 
roll-up and the decay process, cannot be changed independently. For 
example, a change in tip shape changes the core size, as well as the velocity 
and turbulence distributions. High levels of turbulence result in an increased 
diffusion of vorticity, which in turn increases the core size. 


Probe sensitivity of the vortices: Flow visualization studies suggest that 
trailing vortices are extremely sensitive to disturbances created by even very 
small ‘probes or bubbles. This forces one to use non-intrusive means 0 
measurements such as an LDY. Even then, ‘vortex wandering’ [Ref. 8], which 
makes the vortices appear larger than normal in time-averaged velocity 
measurements (for. vortices generated by a wing in a wind tunnel), or the 
unsteady nature of the flow (for vortices generated by a Wing in a tow basin) 
makes the mean velocity profiles in the vortices difficult to determine. 


Large-scale. instabilities: The vortices are seldom observed to decay away 


- owing to viscous and turbulent dissipation, but are, almost always. destroyed 


by either mutual induction instability (Crow instability (Ref. 9]) and/or vortex 
breakdown. The Crow instability grows exponentiallv, and results either ina 
linking of the vortex pair into,a series of crude vortex rings or in a highly 
disorganized intermingling of the. vortices. Vortex breakdown. whose 
Mmeriematcalscefails have mot yet beemeadequately treated. rearranges the 
Vortex structure and increases thé core size. turoulence. and energy dissipation. 
Thus, it is verv difficult to measure accurately the trajectories of the three- 
Gime lonierortices ifolm their creation to their ultimate demise. 


Revnoids number: Even the highest Revnolds numbers, based on wing chord. 
Feaemed i wand tunmels or tovime basins. are in order_of magnatude lower 
than what is possible for an aircraft. Thus, the scale effects are not easv to 
ASSESS. 


Ambient conditions such as turbulence and stratification play major roles in 


the evolution of vortices. The quantification of these eflects requires 
numerical analvsis and extremely careful experiments. 
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(f) Ground or free-suriacevciiecrs. saa. vortex. pair mav move towards a mngid 
boundarv at which.the no-slip condition must be satisfied or towards a free 
surface at which the zero-shear condition must be satisfied. In either case, the 
vortices come under the influence of their images and move accordingly. 

The phenomenon is further complicated by several additional facts. When the 
vortices are propelled towards a ngid surface, vorticity of opposite sign is generated on 
the no-slip boundary and swept towards the vortex pair. The total vorticitv diminishes 
very quickly as vorticity from the two regions diffuses, the wall region serving as a 
strong sink for the vorticity associated with the ‘original vortex [Ref. 10]. The 
developement of a boundary layer along the rigid wall mav give rise to flow separation 
for sufficiently high Reynolds numbers. With or without such a separation, however, 
the center of the vortex pair eventually moves away, or’ ‘rebounds’, from the wall 
(Refs. 10-12]. 

For the case of a zero-stress boundary, the free surface still acts as a vorticity 
sink, but this is relatively weak due to the absence of intense oppositely-signed 
vorticity. Thus, in the absence of other impeding phenomena, one expects a muld 
interaction between the vortices and the free surface and a small rebound of the vortex 
center from the free surface. Howeéver,the ability of the free surface to deforma 
the influence of strain fields leads to a strong interaction between the vortices and the 
free Suriace: | 

It is evident from the foregoing that the motion and the life-span of trailing 
vortices are goverened by a number of nonlinearly-dependent complex phenomena. A 
number of experimental and analytical studies have been carried out at the Naval | 
Postgraduate School by Sarpkaya and his students [Refs. 6,13-17] in ondeniige 
investigate the effects of these parameters on the rise and demise of the trailing vortices 
in homogeneous and density stratified media. These studies have clearly identified the 
various demise mechanisms in both media and established basic relationships between 
the rise of vortices and the governing parameters in a finite as Well as effectively infinite 
medium [Ref. 18], free from ambient turbulence. The present investigation is a 
continuation of the foregoing studies and is confined to the computer modeling. via a 
finite differencing technique, of the rise and demise of trailing vortices in stratified 


media. 


ee Biv AIPEOUIPMIENT AND°PROCEDURES 


a FL OWIRMENT 

The equipment used to generate the trailing vortices has been extensively used at 
this facility over the past three vears [Refs. 6,14-19]. Only the salient features, most 
recent modifications and the adaptation for this work are briefly described in the 
following. | 

The system consists essentially of a towing basin. The auxiliary components 
consist of the brine fill system, turbulence management system, top and bottom 
carriages, velocity measuring system, lighting system, and the model (Refs. 14-19]. The 
brine full system is controlled by a computer to provide a consistently repeatable 
density stratified media [Ref. 19]. . 

Drains are provided at the bottom of each end of the basin. Two parallel rails 
are mounted along the bottom of the tank. A carriage rides smoothly on these rails 
and provides the test body with a constant velocity through the use of an endless cable 
and a DC motor. The velocity of the model is measured and continuously monitored 
through’ the use of a magnetic linear displacement transducer. It yields a signal 
proportional to the velocity of the model with an error of less than one percent. 

The two rails, the carriage and the filling pipes are located on or near the bottom 
of the basin and under a turbulence management system. This system consists of one 


inch thick polyurethane foam sandwiched between two perforated aluminum plates. 


Be MODEL 

A rectangular foil (NACA 0012) was used in the present study. The foil had the 
following dimensions: B = 4.50 in., Ie. = 3.53 ieecnd CC =) 2.32 in, he interior of 
this mcedel was hollowed and used as dve reservoir to seed the vortex cores. The model 
Was mounted on its base by means of a thin streamlined aluminum bar with a cross 
section of a NACA 0006 foil and set at the desired angle of attack. As noted earlier. 
the model was pulled by means of a DC motor, pulley, and cable svstem at the desired 
speed. Additional details of the model construction and mounting are discussed by 
vommson (Ref, 15]. 


bs 


C. TESURrPROGCEDU Ras 

A model is placed in the basin and the basin is filled gradually with a brine 
solution [Ref. 19] to the desired level. In the cases where dye was used to seed the 
vortex pair, the model was filled with dve prior to filling the basin. After removal of 
trapped air and after sufficient time to allow for the escape of dissolved air in the water 
and the elinunation of anv internal currents in the basin, the model was set in motion. 
Most of the experiments were repeated at least three times. 

The motion of the trailing vortices were recorded on highspeed Polaroid film at 
theses. section (ome saleuac plexiglass panels near the middle of the basin). Each 
picture included two clocks accurate to 0.1 seconds, the vertical and horizontal scales 
on the plexiglass window and, of course, the side view of the trailing vortices as they 
rose from the model after formation. The time interval between successive pictures 1s 
determined from the two clocks. The first picture always captured the model the 
instant it passed through the test section. The Subsequent pictures (taken at about 
0.75 second intervals) captured the mse and demuse of the vortices. The vertical rise of 
the vortices is determined from the verticle scale. Attention has been paid to the fact 
that the vortices are farther away from the scale on the window and that the scale 
placed vertically in the middle of the test section does not exactly correspond to the 
scale marked on the window due to refraction and parallax. The necessary correction 
was made by photographing a scale placed in water in the middle of the test section 
together with a scale marked on the window. This resulted in a simple conversion 
table which enabled the determination of the actual position of the vortex from the 
scale reading on the photograph. | 

The results are normalized and plotted in various forms and compared with those 
obtained in the previous runs. Each experiment was repeated three times for most of 
the data presented herein. All trailing vortices were recorded on film until the time 
they have completely dissipated either due to aging (diffusion of vorticitv due to 
viscousity, turbulence. entrainment, and detrainment), or due to Crow instability 
(sinusoidal instabilitv and linking of vortices), and, or due to vortex breakdown (core 
bursting). It was thus’ possible to determine the life spam of vortex cores {rom tiiam 


formation to their final demise. 
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Dew CONMPPUTER MODEL 

The computer model (Appendix A) is run on the IBM 3300 main-frame computer 
at the Naval Postgraduate School. The computer code was written to run on either the | 
MVS Batch Processing system or the VM;CMS svstem. The stratification parameter 
Speeimed [or the expenmiental runs iseused to determune other parameters [Ref. 17] used 
by the program. The code is compiled under VS Fortran and is loaded and executed 
using DISSPLA version 10.0. The code was written such that the user has the capacity 
to Seeulate the mimGer of time step iterations that are executed. Also. the time steps 
at which a constant density contour, complex velocity profile, density perturbation 
contour and the vorticity contour are plotted, are designated by the user prior to 
compilation. These plots are analized to deterrmne the normalized nse height of the 
initial vortex pair until its demise. These are then compared with those obtained 
experimentally. 

The computer model is an intense computational algorithm that requires 
approxiamately 12.5 minutes of Central Processing Unit (CPU) time per iteration. Due 
to this high use of CPU time, the code has been designed to be executed until stoped at 
a user defined step iteration and then save all necessary data required to restart the 
model. The code can then be restarted at that step iteration at a later time. The 
program VORRS (VORtex ReStart) has been specifically written to accomplish this 
task. . 
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Ii. ANALYSIS 


4. NUMERICAL ANALYSIS OF THE VORTE Soi ens 

The generation of the internal waves and the rise and demise of the vortices in a 
stratified medium mav well be analyzed through the use of the equations of motion for 
an incompressible fluid for both laminar and turbulent motions provided that a suitable 
turbulence closure model is adopted and the usual Boussinesq approximation 
gravitational acceleration is much larger than the fluid accelerations) is made. For the 
type of motions considered herein the Boussinesq approximation is quite valid and has 
been used in the investigation of all tvpes of internal waves in stratified fluids. 


For a two-dimensional flow, with y vertical and x horizontal, the equations of 


motion are 
du du du 1 dp 5 
— yy ee eae ee (3.1) 
ot Ox Ov 9 OX 

and 
OV OV OV 1 dp 5 ; | 
— + $ VS = See | (3.2) 
ot OX Oy p dv i 


and the equation of continulty 


Bul oman 
- 


ca = 0. BB 
Ox OV oa 





Defining vorticity as usual by 





and eliminating pressure between Eq. 3.1] and Eq. 3.2 , one has 
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SS et ps £2 ; (375) 
ot Ox OV an erox 


in which p’ is the fluctuating part of the density p given by 
Ome) ap (X.Y,0) (3.6) 


Where pis the reference density, and p(y) the initial density at elevation y att = 0. 
The diffusion of density is given by 
0 0 CG : 
+us-—+v— — vV*)p = 0 (oa7) 


orn OX Oy 


Combining Eqs. 3.3. 3.6 and 3.7, and simplifving, one has 


Cl 


dup’ dvp’ g Gs 
my oP, oe ane vv-p) + v— 
Gt OX Oy Ov OV 


a 


(3.8) 





ty 


It is convenient to cast the foregoing equations into nondimensional forms, 
scaling each variable by a quantity characteristic of its expected magnitude. Two 
massitele time scales exist. There is the dynamic time scale, which is the time a 
characteristic length would be traversed by a fluid particle traveling at the characteristic 
velocity, and there is the buoyant time scale based on the natural buoyancy frequency 
of the stratified flow, 1.e., the Brunt-Vaisala frequency N, as defined by Daly [Ref. 17]. 
Each of the scales gives a slightly different form of the normalized qoverning equations. 

Dynamic Scale: | 3 
Introducing U. and L. as the characteristic velocity and length, one has 








ei Lea ae L 
¢ = S to=-* Re=—S F = —<, 3.9) 
a aie i v8 Lg 7 





at Gu. me lL ¢ 
mM + (em + mm) = __ GEG fp (3.10) 
at, Ox, OY Re ee , 


and 




















Cn’ Gu  6v 670 
on Ox. Oe Re OY 
earahy a yietal 

: ie Uae < 
F- = — Re = —— +, n = — 

. oL V N. 

Saw C os 

w2e% y2-- SPM _ 8B _w 

: e : p, dy oan e. 


~ 
A 


Equations 3.10 ans 3.11 are valigewiener = tae when the bouvancy has 
very little influence on the nonlinear dynamics of the motion. In this case, the density 
perturbation acts as a passive scalar advected by the velocity field. 

Bouvant scale: 


Introducing the following dimensionless parameters: 


= p u V X 
Pin —— = u SS ——, we = eae >< Se 

Py ™ U, i U, ™ L, | (3.12) 
Mpa a {=e ca = = N? = 

Ls ‘, L. 


and normalizing Eqs. 3.5 and 3.8, one has 


oC du_¢ vam le , 

=irinee F (—2 HG et mom ) =- _¥G4 aa oP (3. 

ot, OX One Re Geos 

and 

op du_p Gv_p 

—M + P(_mem + mm) = F p*y 3. 

a ex ay en - 
mM ae ¢ 0 | 


in which F,, Re and n° are as defined in Gen oaiele 
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Equations 3.13 and 3.14 are valid when F. < < 1 and buoyancy dominates the 
flow. When F. approaches zero, the equations that result from Eqs-3.13 and 3.14 
describe the propogation of linear internal waves. The F, << 1 regime is of interest 
in the present investigation because for submerged bodies of naval interest F. (the 
Froude number) is about 0.01. The analvsis will consider the full nonlinear equations 
Biempoy Eds. 5))> ard 5-12 rather than their linearized form (1... the case of F. = 0). 
However. the flow will be assumed to be inviscid, 1e., the viscous diffusion will be 
ignored to a first order approximation. Subsequent analysis will incorporate the effects 
of viscous and turbulent diffusion into the numerical calculations. The velocities u and 


v are given by the Biot-Savart law as 


(y’ - ¥)G(x’y )dx’dy’ 


u(X,Vv) = Sale 
(X,y) J a (oat) 
and 
(x = X)E(X’ y')dx’dy’ 
———— (3.16) 


in which (x,y) is the point at which the velocity is calculated and (x’,y’) is an arbitrary 
point at which the vorticity is ¢(x’,y’). 
The radial distance r is given by: 


i (X - oe + (yy - Pe (3.17) 


Normalizing Eqs. 3.15 - 3.17 through the use of the characteristic dimensions 


eiren in Eq. 3.12 


Dee eee ex dy , : 
eu. - 7 (3.18) 
— poi oe Ge Sm ‘dx dv 3.19) 
en on? oli 


and 


po? S(xe =k) Ge | (3.20) 


It is appropriate to take V., the initial mutual induction velocity of the vortex 
pair, as the characteristic velocity and the b,, the initial separation of the vorticies. as 


the characteristic length. Thus. one has: 








L.=,. U,=¥, =F = (3.21) 
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in which aoe is called the stratification parameter SP. The value of SP normally 
varies from zero to about 100. The larger the SP, the stronger the stratification. For 
example, for a model running at a speed of 5 ft'sec at an angle of attack (a) of 10°, a 
SP smaller than 1.0 represents weak stratification, a SP between | and 5 represents 


medium stratification, and a SP larger than 5 corresponds to strong stratification. 


B. COMPUTER AIDED MODELING OF VORTEK MOTION 

The computer algorithm for the numerical solution of Eqs. 3.13 and 3.14 is 
provided in Appendix A. The solution is done assuming the case of inviscid flow and 
the initial vorticity distribution was assumed to be Gaussian, 1.e.: 


< r 2 


, | F : . : - 
Ca y exp(— rie x) (3.23) 
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A computational domain of 5b, by 10b, is used in the algorithm in order to 
mumick the physical features of the tow basin. The no-slip condition is satisfied an the 
walls corresponding to the chanel walls. The zero normal velocity condition is assumed 
at the free surface and along the vertical line of flow symmetry. All computations are 
done in the fourth quadrant. The fields of velocity, vorticity, density, argues 


fluctuating components of the density are calculated at each time step and plotted at 


regular intervals. The computer algorithm is quite robust and can be applied to any 
stratified and unstratified flow situation. As the velocity field is computed from 
analytic expressions (i.e., the: Biot-Savart equations) vice an algorithm involving an 
iterative scheme, stability of the program is enhanced at a slight expense in execution 
time. An upwind finite differencing technique [{Ref. 20] was used with an original grid 
dimension of 24 by 44 nodes with a spacing of 0.25b,. This grid size was subsequently 
increased to 34 by 64 and 44 bv 84 in order to exanune the dependency of the results 
@aeene Crid size. 

The mesh coordinates (x,v) covering the computational domain (of dimension 


m-4 by n-4) are: 


ee | 2a daw tiie 2) 
(3.24) 
Y € {2,3,4,...,n-2}. 
The necessity for the the overlap of two nodes on all sides of the computational 
domain is to facilitate satisfving the boundary condiditons. Eqs. 3.13 and 3.14 are 
solved over the original mesh which results in vorticity and density information. 

The two complex velocity grids are derived by shifting the original mesh h/2 in 
tae. xX direction for the u velocity and h/2 in the y direction for the v velocity. The 
values of these new nodes are determined using the vorticity derived’ for the original 
‘computational grid using the Biot-Savart law, Eqs. 3.15 and 3.16. The velocities at the 
shifted grid points are then averaged across the computational grid points for further 


use in the finite differencing scheme. 


IV. DISCUSSION OF RESULTS 


The computer code described in the previous section has been used to carry out 
calculations for four specific cases. The first of these dealt with the case of the rise of 
vortices in a homogeneous unstratified medium. The grid size as well as the core 
radius were varied within limits to explore the differences in the motion of the vortices. 
Figures B.] through B.I3 show the velocity and vorticity contours at various 
normalized times. Clearly, the initial vortex spacing increases with time and the 
vortices remain symmetrical due to the conditions imposed on the numerical analysis. 
In other words, these vortices cannot be subjected to sinusoidal instability or core 
bursting. Thus, the numerical results should agree more closely with the experimental 
cases Where the vortices did not undergo such instabilities. Figure B.14 shows the 
experimental and numerical results. The experimental data are separated as ‘HIGH’ 
and ‘LOW’ because of the fact that the crests and troughs of the sinusoidal instability 
‘Were separately evaluated. The numerical results for three normalized core sizes are 
also shown in Figure B.14. It 1s clear that the core size has a pronounced effect on the 
numerical predictions and that the results obtained with r, = 0.08 agree more closely 
with the experimental data. As noted earlier, the numerical results should be compared 
with the average instantaneous height rather than with the ‘HIGH’ and ‘LOW’ of the 
vortices. There is no conclusive analvsis or data in the literature regarding the vortex 
core radius, It has been shown previously (Ref. 16] that the core radius "@epemee 
strongly on the shape of the trailing edges of the lifting surface. The better rounded 

the edges of the control surface the tighter the winding of the vortex sheets and thus 
the smaller the vortex core. It appears that a vortex core of 0.07 or 0.08 should be 
used in future calculations. | 

Figures B.I5 through B.351 show lines of constant density, velocity, density 
perturbations, and vorticity contours for the statification parameter of 1.0. FL = 
0.02976, r, = 0.09. Y’b, = 78.0vand X-b, =O: The gmd used wase4siieaie 
These figures show that at small times the vortices rise vertically and the circulation in 
the flow field 1s primarily due to the initial vortices. As time increases {see €.9))tamme 
B.21b). the regions of circulation with countersigned vorticity develope in thewiigen 


right and lower left regions of the vortex. It is easy to show that this counter vorticity 


not only reduces the rise velocity of the original vortex pair but pushed the pair against 
each other. Consequently, the countersigned vorticity begins to dominate the flow (see 

g., Figure B.28b) and the vortex pair mugration stops. With further increases in time 
the vortex pair begins to nugrate downward. 

The density contours reveal the same phenomienon in a different context. As the 
emmices tise, fluid Of Greater density is pushed upwards (see e.g., Figure B.25a) into 
regions of lesser density. Since such a migration cannot go on indefinitely, the vortices 
rise tO a maximum height and then begin to sink downwards. The calculations do not 
take into account thé sinusoidal instability and the vortex breakdown. Consequently, 
the vortices in the numerical calculations continue to exist unimpeded by the instability 
mechanisms. In realitv, the vortices begin to break up as they near the end of their 
maximum migration and eventually disappear. 

The experimental and calculated values are compared in Figure B.32. The 
correspondence between the measured and calculated values is surprisingly good up to 
the time of maximum rise. This is partly because of the experimental fact that the rise 
of vortices in a highly stratified medium ee OO) = 0.02976) is not 
strongly dominated by sinusoidal instability or vortex breakdown. The migration of 
vortices is inhibited primarily due to the reduction of vorticity of the initial vortices and 
the creation of countersigned vorticity. 

The numerical calculations were also applied to a medium with a relatively 
smaller stratification for which the rise of vortices are closer .to that of the 
homogeneous case and thus the effect of the core radius is more pronounced. The 
results obtained with ‘oe = 0.50, P, = 0:02976 and with es 0.08 are shown tn 
Figures B.33 through B.43. The observations made earlier regarding the velocity, 
vorticity, and density contours may be made with this particular stratification also. 
Obviously. the vortices rise to a greater height and the countersigned vorticity is 
relauvely weaker and developes at a later stage. The calculations with this particular 
stratification have been repeated by varving only the initial core radius. Figure B.44 
shows the representative experimental data and the numerical results obtained with 
initial core-radi of 0.07, 0.08. and 0.09. Evidently and as anticipated on the basis of 
the previous calculations the initial core radius of 0.07 vields results which are in 
excellent agreement with those obtained experimentalv. Figure B.44 also shows that 
the rise and denuse of vortices must be a strong function of their initial core radii. 


Experiments in both the laboratory and the field have shown that to be true and 


pointed out the fact that the life of a vortex is dictated partly by the conditions 
prevailing at its creation and partly by the conditions surrounding the vortex at later 
stages of its motion. | 

The results presented so far have dealt with vortex pairs released at relatively 
large depths of submergence and have shown that the computer code predicts results 
which are in good agreement with those obtained experimentally provided that a core 
radius of about r, = 0.07 is used in the calculations. It has further been shown that 
the smaller the stratification parameter the more dependent the results become on the 
initial core readius. Thus, the numerical analysis has not only provided a power of 
prediction for the motion of the vortices but also helped to delineate the proper value 
of the imitial core radius. 

Having assertained the characteristics of the model and the features of the 
velocity, vorticity, and density variations it was deemed necessary to explore problems 
of greater practicle and operational concern. It is because of this reason that the depth 
at which the vortices are generated was chosen to be Y/bo = — 2.0 and the 
subsequent motion of the vortices was explored. Figures B.45 through B.52 show the 
density, velocity, density fluctuations and vorticity contours at suitable time intervals. 
As the vortices approach the free surface and begin to move sideways under the 
influence of their images, the constant density contours begin to intersect with the free 
surface (see e.g., Figure B.49a) and move sideways also. The intersection of these 
density contours is expected to give rise to internal-wave-generated surface scars. The 
prediction of the characteristics of-these scars in terms of the vortices is of importance 
and requires that the free surface condition be changed. This change is to allow for 
free surface deformation under the influence of the tangential as well as normal 
velocities induced at the free suridee oy tieeveruces, 

- Additional calculations are underway to examine the effects of free surface 
deformation, viscous diffusion and ambient turbulence for various degrees of initial 
Stratification. The results presented herein and the initial comparisons are extremely 
encouraging and are expected to lead to a better understanding of this eXtremiam 


complex and challenging phenomenon. 


Pee OC CUSIONS 


The investigation reported herein warranted the following conclusions: 
(1) The motion of trailing vortices generated bv submerged bodies in a stratified 
Or unstratified medium can be effectively modeled bv the developed computer 


algorithm: 


(2) The computer model has aided in the delineation of the proper value for the 


initial vortex core radius; 


(3) The initial vortex core radius and vorticity strength are the primarv factors 


governing the stability of the numerical results; 


(4) The migration of vortices is inhibited primarily due to the reduction of the 


vorticity of the initial vortices and the creation of countersigned vorticity; 


(5) The computer model could be made more realistic through the inclusion of a 


deformable free surface in the computations. 
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APPENDIX A 
VORTEX WIOTION CODMIPU TE ate err 
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PROGRAM VORTEX 


ZQ NOVEMBER 86 PROGRAMMER : BREAN S.eLLEER 


ASSAY) “On “Oo 


CRAKKAARKAKAKKRKKARARRERRRERRRRKRKRRRAKRRREKRKRRRRKRKRKRREKRRKRKRERKRRRKRKRKRARRKARR 


C 

C PURPOSE: 

C 

C THIS PROGRAM NUMERICALLY SIMULATES THE RISE OF A TRAILING 

C VORTEX PAIR SHED OFF A LIFTING SURFACE SUBMERGED IN A 

C DENSITY STRATIFIED MEDIUM. 

i  s— 
C 

C TO RUN THIS PROGRAM: 

C 

C 1) GO TO "INPUT PARAMETERS" SECTION OF CODE AND CHANGE ANY PROBLEM 
C PARAMETERS ( SP, FV, DT, ZTMAX, RO, MP, NP, ISAVE) AS DESIRED. 
© } . | 

C  .2) SET MAXIMUM TIME IN DATA STATEMENT MAXSTP; 

C EX: DATA MAXSTP /N/ 

C ( WHERE N= MAXIMUM NUMBER OF TIME STEPS) 

C 

C 3) SET INTERVALS FOR PLOTTING BY ENTRIES IN DATA STATEMENT NPLOT; 
C EX: DATA NPLOT/ 1,2,8,12,20,40,50,75/ 

C ( THIS PLOTS OUTPUT AT THE 1ST,2ND,8TH,12TH, 40TH, 50TH, AND 
C 75TH PROGRAM TIME STEP) 

C 

C 4) CHOOSE WHAT OUTPUT DEVICE TO USE, EITHER TEK 618 OR COMPRES 
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(METAFILE OUTPUT) BY GOING TO GRAPHIC DEVICE SECTION OF PROGRAM 
AND COMMENTING OUT WHICH OF TWO DEVICES YOU DONT WANT TO USE. 


5) COMPILE PROGRAM UNDER VS FORTRAN USING THIS COMMAND; 
BOR LV Ss a) ORE x (OPT(3) 


(Noite THise OPTIMIZES GQBE AD LEVEL 3 ) 


6) GO TO APPROPRIATE TERMINAL. .AND RUN UNDER DISSPLa. 

Pee nee Disor hae 2AEC PROMPTS YOU ASK FOR 5S CYLINDERS OF 
Derren ane DISK STORAGE ANDWAESO YOU MAY WISH TO ISSUE 

5B. (WHEE peees EC BPROMPERS YOUNYOUNMAY WISH TO ISSUE: 

FILEDEF 06 DISK VORILOUT LISTING (PERM 

Ce none newnoe taee ra tira OUTPUT TO A DISK FILE 
NAMED: VORLOUT LISTING Al 

Cee ee eel ORNS ENDING OUTPUT TO YOUR DESIGNATED 
DERIIR@ES.. 


ADDITIONAL PROGRAM NOTES: 


er) Tails PROGRAM COMPUTES DENSITY AND VORTICITY IN A (M-2 K N-2) GRID 
Bt SOLVINGPTEE BOUYANETLY S@=LED FORM, OF THE NAVPER-STOKES 
POUR GNS Ven AN “UPWOIND DIFFERENCING"’ FINITE DIFFERENCE METHOD. 


2) THE VELOCITY FIELD IS COMPUTED USING A FORM OF THE BIOT-SAVART 
LAW AT ZACH NODAL POINT INCLUDING THE CONTRIBUTION OF THE 
So ctevOnneeitS “ADDED TO INCLUDEATHE FREESSURPACE EFFECT. 


3) THIS VERSION OF THE PROGRAM INCLUDES A DISSPLA ROUTINE 
TO PLOY THe DENSITY “PERBUBATLONS=, CONSTANT DENSITY» 
eet ie nneerry FLED AT EACH SELBG@TSD TIME STEP. 
Tipo SekLeChiONSMEORePLOTTONG ARE INDICATED BY “EEMERIES 
IN THE ARRAY NSTEP. MAXIMUM NUMBER OF TIME STEPS IS INDICATED 
BY AN ENTRY IN THE ARRAY MAXSTEP. 


Lo 
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4) THE COMPUTATIONAL GRID (M-2)X(N-2) IS USED FOR ALL COMPUTATIONS - 
EXCEPT FOR VELOCITY CALLS TO THE BIOT-SAVART SUBROUTINES. 
TWO ADDITIONAL GRIDS, ONE FOR U AND ONE FOR V, ARE SHIFTED #/2 
UNITS IN THE X AND Y DIRECTION RESPECTIVELY. THE SUBROUTINES 
COMPUTE U AND V VELOCITIES ON THESE GRIDS USING VORTICITY FROM 
COMPUTATIONAL GRID. THE VELOCITIES ARE THEN AVERAGED ACROSS THE 
COMPUTATIONAL GRID POINTS TO BE USED FOR THE FINITE DIF FERENGMMG: 


5) THIS VERSION OF THIS PROGRAM WILL BE WRITTEN IN, SINGLE PRECISION 
AND IN AN OPTIMIZED FORM OF THE ORIGINAL CODE. 


6) THIS PROGRAM HAS THE OPTION OF SAVING THE REO AND ZETA ARRAYS 
AND OTHER PROGRAM CONSTANTS IN ORDER TO RESTART THIS SIMULATION 
AT THE TIME A PREVIOUS RUN HAD FINISHED. PROGRAM VORRS [5 
SPECIFICALLY WRITTEN TO INPUT THE DATARPTSEee ete BY THis 
PROGRAM AND RESTART THE SIMULATION WITH THE SAME PARAMETERS. 


TO ACCESS THIS SAVE FEATURE SET VARIABLE ISAVE =1 . A FILEDEF 
OF FORM ... FIGEDEF 08 DISK <BND> <FT> <ai> ... MUST BE iseame 
AT THE APPROPRIATS Cite” IN Eemeon ie. 

(NOTE: TO DISABLE THIS SAVE FEATURE SET ISAVE =0 ) 


7) TO GHaNGs THE MESH"DENSITY OF THE COMPUTATIONAL. DOMAPR [He 
USER SHOULD D@P THE FOLLOW Itic : 
A. ADJUST M AND N 
ADJUST INITIAL VORTEX POSITION: Die’, “eee 

C. ADJUST THE DIMENSIONS OF THE WN ARRAY IN SUBRAUeE TE: 
CPLOT, ZPLOT, AND DRPLOT TO BE (MW,NW) 

D. ILE PRE SGRED Tseb33ce THAN 30 X 100 THEN ADJUST THE 
DIMENSION STATEMENTS IN THE MAIN PROGRAM AND IN 
SUBROUTINES ZPL, VPLOT END VZrn0l TOS trae ae es 
DIMENSION. 


e 
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MAJOR VARIABLES USED IN THIS PROGRAM; 


Poe VARTSSLES; 


cy = NONDIMENSIONAL TIME 

DT See SLEE 

H = NONDIMENSIONAL GRID LENGTH PARAMETER= DX = DY 
RO = NONDIMENSIONAL VORTEX CORE RADIUS 


XWIDTH= NONDIMENSIONAL WIDTH OF COMPUTATIONAL AREA 
ZTMAX = MAXIMUM VORTEX STRENGTH 


Se = NONDIMENSIONAL STRATIFICATION PARAMETER 
FV = NONDIMENSIONAL FROUDE NUMBER 
BVND = NONDIMENSIONAL BRUNT-VALSA FREQUENCY 


INTEGER VARIABLES; 


M = NUMBER OF NODES IN X DIRECTION 

N = NUMBER OF NODES IN Y DIRECTION | 

MP = X NODAL POINT OF VORTEX START POSITION 
NP = Y NODAL POINT OF VORTEX START POSITION. 
MM = X ARRAY DIMENSION 

NN = ¥ ARRAY DIMENSION 


NFIG = PLOT NUMBER 
ITERATION NUMBER 


ra 
W") 
4 
mH 
‘d 

ll 


ARRAY VARIABLES ; 

RHO = DENSITY PERTUBATIONS AT EACH NODE DUE TO DYNAMICS 
RHOT = INITAL VALUES OF DENSITY AT EACH NODE DUE TO GRADIENT 
Pea — cUMeOrSSOTH DENSDIY ERFECTS AT EACH NODE 

RENEW = NEW VALUE OF RHO AT NEAT TIME STEP 

Zoe Venn tL CITY Qe ACH MODE 
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CTE = ite VALUE OF VORTICITY AT MEXT TIME STEP 


Cyr ee), 3G “Gy C2 Sere) ies FO). Ee 


M 


= X COMPONENT OF VELOCITY AT EACH NODAL POINT 
v = Y COMPONENT OF VELOCITY AT EACH NODAL POINT 
U2 = U VELOCITYS COMPUTED AT BIOT-SAVART NODAL POINTS 
V2 = V VELOCITYS COMPUTED AT BIOT-SAVART NODAL POINTS 
= X COORDINATES AT EACH NODAL POINT 
= Y COORDINATES AT EACH NODAL POINT 
XU «= KX COORDINATE AT EACH U 8. S$. NODAL POINT 
YU = ¥ COORDINATE AT ZACH U B. S. NODAL POINT 
XV «= -X COORDINATE AT EACH V B. S. NODAL POINT 
YV = Y COORDINATE AT EACH V B. S. NODAL POINT 


SUBROUTINGS US2eD Airis heen ie 


SUBROUTINE COUT OUTPUTS NODE COORDINATES 


SUBROUTINE OUTPT 


OUTPUTS VARIOUS PROGRAM ARRAYS WHEN CALLED 


SUBROUTINE CPLOT CREATES ARRAY OF DENSITY AT EACH POINT 


AND CALLS DISSPLA CONTOURING ROUTINE TO PLOT. 


CREATES ARRAY OF VORTICITY AT EACH POINT Same 
CALLS DISSPLA COUTOURING ROUTINE [TO Pitas 


SUBROUTINES Zr Eer 


SUBROUTINE DRPLOT- CREATES ARRAY OF DENSITY PERTUBATIONS AT EACH 
POINT AND CALLS DISSPLA CONTOURING ROUTINE 
TOMELOm: : 


SUBRGURE ize rt - OUTPUTS VORTICITY CONTOURS AS A PRINTER tee 


SUBROUTINE VPLOT - CREATES COMPLEX VELOCITY VECTORS AT EACH 
NODAL POINT AND THEN CALLS THE FOLLOWING 
SUBROUTINES TO PLOT THEM ON DISSPLA: 
A) V2PLOT 


ee ree Gn SC en (eG 0) «69°C 0. 8 OY I I OO a) 6 


Cle leet eee) (2 Cd ete (Cy) (2. ( CY CY OO OO OC 


B) AMIN 
e) AMAX 


See ee meme Loree a LTO CONTOUR VORTICITY ON PRINTER. 


PONG@RPONSSUSED IN THIS PROGRAM: 


etc Oem mecOUrUGES THe VY COMPONENT OF VELOCITY USING 
THE BIOT-SAVART LAW. 


MeertON Uss - COMPUTES THE U COMPONENT OF VELOCITY USING 
Ree BLOL= SaVART Lay. 


ere i lemmr NEW = TAKES THE FINITE DIHBERENCE TO STEP THE 
EQUATIONS IN TIME. 


“DIMENSION setisc, 100) ,U(50,100) ,V(50,100) ,ZETA(50,100) 
,ZTNEW(S50,100) , RHNEW(50,100),NPLOT(8) ,X(50),¥(100) ,RHO(50,100) 
,RU(50), ¥U(100), xV(50), ¥V (100) ,U2(50,100), v2(50,100) 

CHARACTER*8 LABEL 

COMMON WK(15000) 

COMMON /PN/ DT,1,J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT, FV 

COMMON /PARM/ H;M,N,PI 

DATA MAXSTP/40 / 

DATA NPLOT/1,10,20,30,40,50,60,70/ 


Gael COMPAS 
GALL TEK612 


MM= 50 
NN= 100 


oP) 
oP) 


DT= 2.000E0 
T=0.0 

XWIDTH = 5.0 

RO= 0.09E0 

SP= 1.20 

FV=0.029765 

ZTMIN= 1.0E-8 

M= 44 

N= 84 

ISAVE= 1 
ee: DESIGNATE NODAL POINTS TO PLACE INITAL VORTEX AT (MP,NP) 

MP=6 

NP=66 


------ COMPUTE ADDITIONAL PARAMETERS --------------------------------- 


..... DENSITY GRADIENT IN Y DIRECTION INITALLY = DRHO 
DRHOI=-SP*SP*EV*EV 


Fae NONDIMENSIONAL B. V. FREQUENCY = BVND 
BVND= SP*SP*FEVAEFV 
are MAXIMUM VORTICITY AT ANY POINT = 2TMaAx 


ZTMAX=-(FV/(RO*RO)) 
. GEOMETRIC CONSTANT PI 
PI= 4.0E0*ATAN(1.0E0) 
ea. GRID LENGTH PARAMETER = H 
H= XWIDTH/(M-4) 


Sgt COMPUTE MISC M AND N PARAMETERS 
Mee= M=1 
MES hal 
Mit2= M2 
NP1=N+1 | 
NM1=N-1 


NM2= N-2 
MW = (2*M) - 7 
NW = NM2 - 1 


onl pee ee OO oes more TGINS “<< - omen renee nnn nr ernrren= 


X(1)= -H 

Y(l)= 4H 

XU(1)= -1.5E0*H 
YU(1)= H 
XV(1)= -H 


eyGy= 1.5E0*H 


-------- COMPUTE COORDINATES ------------------ 
( NOTE: COMPUTATIONAL GRID IS IN 4TH QUADRANT +X -Y ) 


BOmto. I= 2,MPi 
— XU(I)= XU(I-1)+H 
XV(I)=XV(I-1)+ H 
10 Mé)o= X(P=1) + H 
DO 11 J = 2,NP1 | 
YU(J) = YU(J-1) - HD 
YV(J) = ¥V(J-1) - H 
a way =" Y( 3-1) - H 
| YHT = ABS(¥(N-2)) | a 
C...... OUTPUT INITAL PARAMETERS OF THIS RUN OF PROGRAM.....ccsececeees 


WRITE (6,1275): 
WRITE(6,1500) 
WRITE(6,1535) DT,H 
WRITE(6,1540) XWIDTH,YHT 
WRITE(6,1550) SP,FV 
WRITE(6,1560) M,N 
WRITE(6,1570) BVND 
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(ae OUTPUT COORDINATES OF NODES..... (os sade dee 
WRITE(6,1580) X(MP),Y(NP) 
WRITE(6,1590) 
CALL COUT(X,Y,MM,NN) 


( FOR GRID SYSTEM X,Y COMPUTATIONAL GRID IN THE 4TH QUADRANT ) 
( TERMS ARG1 - ARG4 ARE IN QUADRANTS 1 - 4 ) 


ARGMN= 30.0 
XCENT= X(MP) 
YCENT= Y(NP) 
DEN= 2.0E0*RO*RO 
ZETA(MP,NP)= ZTMAX 
DO 30 J=1,N 
DO 30 I=1,M 
TRM1=0.0E0 
'TRM2=0 . 0EO 
TRM3=0 .0EO 
TRM4=0 .0E0 | 
ARG1= ((X(I)-XCENT)**2 + (¥(J)+YCENT)**2)/DEN 
IF (ARGL .GT. ARGMN) GO TO 12 a 
_TRM1= EXP(-ARG1) 
Wy ARG2= ((X(I)+XCENT)**2 + (¥Y(J)+YCENT)**2)/DEN 
IF(ARG2 .GT. ARGMN) GO TO 14 
TRM2= EXP(-ARG2) 
14 ARG3= ((X(1I)+XCENT)**2 + (Y¥(J)-YCENT) **2)/DEN 
IF(ARG3 .GT. ARGMN) GO TO 16 
TRM3= EXP(-ARG3) 
16 ARG4= ((X(I)-XCENT)**2 + (Y(J)-YCENT)**2) /DEN 
IF(ARG4 .GT. ARGMN) GO TO 18 
TRM4= EXP(-ARG4) 
18 IF((J.EQ.NP) .AND. (I .EQ.MP)) GO TO 20 
ZETA(I,J)=ZTMAX*(-TRM1 +TRM2 -TRM3 +TRIM4) 
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Pai Satgeeted)) Lt. ZIMIN) ZETA(I,J)=0.0E0 


20 RHOI(I,J)= DRHOI*Y¥(J) 
RHO(I,J)= 0.0E0 
30 CONTINUE 
WRITES (6,1250) 
C CALL OUTPT(ZETA,'ZTA(I,J)!,MM,NN,O) 
@ CALL OUTPT(RHOI,'ROI(I,J)!,MM,NN,O) 
C 
2 INDDABIZE TIME STEP COUNTERS...... 
NFIG=0 
NSTEP=0 
C 
CRARKRAKAKKRR ENTRY POINT OF TIME LOOP KRRARRKRKKAKARR 
50 CONTINUE 
2 
¢,  _——e COMPUTE U VELOCITY FOR ALL NODES IN U VELOCITY GRID..... 
E 
po 70 I= 1,MP1 
po 70 J= 1,N 
70 U2(I,J)= UBS(XU(I),¥U(J) ,MM,NN,X,Y, ZETA) 
C ) | | | 
C.........-. COMPUTE V VELOCITY FOR ALL NODES IN V VELOCITY GRID..... 
C 


DO 100 I= 1,M 
DO 100 J= 1,NPL | | 
100 V2(I,J)= VBS(XV(I),¥V(J),MM,NN,X,Y, ZETA) 


C 
i ae INTERPOLATE VELOCITIES FROM VELOCITY GRIDS TO NODAL 
C POINTS ON COMPUTATIONAL GRID. 


VMIN= 1.0E-08 
Bemmz0 I= 1,é 
DO 120 J= 1,N 
ie we elie (Il ee Wek I41,35))/2.0£0 
Volo )= (VXI G)\+ Vell, J+1))/2.0£0 
TA w@ASS (WMI ,3)) .LT. WMIN) U(I,J)= 0.0E0 
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IF (ABS(V(I,J)) .LT. VMIN) V(I,J)= 0.0E0 


120 CONTINUE 
IF (NSTEP .£0. 0) CALL OUTPT(V,'V(I,J) 
IF (NSTEP .EQ. 0) CALL OUTPT(U,'U(I,J) 


eo eee#e»eoeeee#eeeeeee#e#ee#see¢¢ ¢ @¢ @¢@6¢@hUcehUhrmhUCUchhUClhhUCUh MhUCUchMhUCUcr HmhUhUhOOCUhHMhUh]]HU UC H HO HF GC 8 ee © ee eee helm 


' MM,NN,1) 
1 MM,NN,1) 


eeee?e?eeeee?e#seee#e# @® @¢ @ @ @ @ #¢ # @ @® @ 


INCREASE Titk COUNTER sv Off@m 
CHECK TTERATION COUNTER AND TERHINSTS LP NSIEBe [5 N@RE Seales 
COMPUTE VORTICITY AND DENSITY AT ALL PNRSRteR rem 
( STORED AS ELEMENTS OF ZTNEW AND RHNEW ) 


NSTEP=NSTEP+1 
IF (NSTEP .GT. MAXSTP) GO TO 1700 
T= T+DT 
DO 600 I=2,MM1 
DO 600 J=2,NM1 
VPT= V(I,J) 
UF= (U(I+1,J)+U(I,J))/2.0E0 
UFA=ABS (UF) 
UB= (U(I-1,J)+U(I,J))/2.0E0 
UBA=ABS (UB) 
VF= (Vil, gel) +v( Gey cee 
_VFA=ABS (VF) 
VB= (V(L,Jt1)+V(I,J))/2.0E0 
VBA=ABS (VB) | 


eeee#eeeeee#e#@ee#?e? 0 08 @¢ ¢ @ @ ¢ @ @¢ @ 


ZTNEW(I,J)= PNEW(ZETA,RHO,MM,NN,-1.0E0,0.0E0,0.0E0) 
600 RHNEW(I,J)= PNEW(RHO,RHO,MM,NN,0.0E0,0.0E0,BVND) 


ee ASSIGN NEW VALUES OF RHO AND ZETA TO 


RHO AND ZETA ARRAYS 


DO 700 I = 2,MM1 
DO 700 J = 2,NMl 
ZETA(I,J)= ZINEW(I,J) 
700 RHO(I,J)=RHNEW(I,J) 
eee, CALCULATE RHO AND ZETA VALUES AT I=M-2 AND J=N-2 BOUNDRIES 


DO 800 I= 2,HM2 

DO 800 J= 2,NH2 

ZETA(MM2,J)= (ZETA(MM1,J)+ ZETA(M-3,J))/2.0£0 

ZETA(I,NM2)= (ZETA(I,NM1)+ ZETA(I,N-3))/2.0E0 

RHO(MM2,J)= (RHO(MM1,J)+ RHO(M-3,J))/2.0E0 

RHO(I,NM2)= (RHO(I,NM1)+ RHO(I,N-3))/2.0E0 
IF(ABS(ZETA(I,J)).LT. ZTMIN) ZETA(I,J)=0.0E0 
IF(ABS(RHO(I,J)).LT. ZTMIN) RHO(I,J)= 0.0E0 

300 CONTINUE 


PLOT FLOW PATTERN WHEN NSTEP = VALUE SPECIFIED IN NPLOT 


THEN RESTART ALGORITHM AT LINE 60 FOR A NEW TIME STEP 
DO 1100 K= 1,8 
IF (NSTEP-NPLOT(K)) 1100, 1200, 1100 
1100 CONTINUE ees as -. : 
GO TO 60 | 
1200 NFIG=NFIG+2 
TM= TREV 


IF (NSTEP .£Q. MAXSTP ) CALL OUTPT(ZETA,'ZTA(I,J)',MM,NN,NSTEP) 
IF (NSTEP .EQ. MAXSTP) CALL OUTPT (RHO, .'RHO(I,J)/,MM,NN,NSTEP) 
WRITE(6,1300)NSTEP,T,TM 

CALL CPLOT(RHO,RHOI,MM,NN,NSTEP ,TM,MW,NW) 

CALL VPLOT(U,V,X,Y¥,HM,NN,TM) 

CALL DRPLOT(RHO,MM,NN,NSTEP,TM,MW,NW) 

CALL ZPL(ZETA,MM,NN,XWIDTH, YHT, ITER) 


SNe 


£256 
yates 
1390 
1400 
1500 
1530 
Le 
1540 


1356 


1560 


S70) 
Te 0) 


1338 


1800 
1900 


CALL ZPLOT(ZETA,MM,NN,TM,MW,NW) 
FORMAT(// 10X,'INITAL ZETA AND RHO VALUES 
FORMAT('1') 


7) 


FORMAT (/10X, ‘ITERATION NUM:' 16,4)" 2iih= eo. ST <= eee 


WRITE (6,1500) 

FORMAT( ///) 

FORMAT(' ',5X,!INITAL PROGRAM PARAMETERS:' ) 
FORMAT('0',5X,'TIME STEP= ',F8.4,5X,'GRID SPACING H = 
FORMAT('O',5X,'NON DIMENSIONAL CELL DIMENSIONS 
* F8.4,5X,'YDIRECTION =',F8.4) 
FORMAT('0',5X,'STRATIFICATION PARAMETER =! 
* 'FROUDE NUMBER =!,F8.4) 
FORMAT('0',5X,'NUMBER OF NODES IN X DIRECTION =',14,5x, 
* ‘IN Y DIRECTION =! ,I4) 

FORMAT('0',5X,'NONDIMENSIONAL B. V. FREQUENCY =',F15.9) 


, Oras Ky 


FORMAT('0',5X,'X COORDINATE OF VORTEX START POSITION =',Fa.4/5m5 


x '¥Y COORDINATE =',F8.4) 
PORMAG (1 0" (ox OuibeUin COORDINATES OF ALL OTHER NODES: ' ) 


LOOP BACK TO STARTING POINT OF ITERATION 
GO TO 60 | | 
CALL DONEPL 
CONTINUE 


ea wes a= = ss . 


IF (ISAVE 
WRITE(38,*) 
WRITE(8,*) 
WRITE(8,*) BVND,H,ZTMAX 
DO 1800 I= 1,M 
DO 1800 J= 1,N 

WRITE(8,*) ZETA(I,J), RHO(I,J) 
STOP 
END 


.EQ. 0 ) GO TO 1900 
SP ipeieve 10 
M,N,MP,NP 
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Peirce) 
; ADTRECT IONS 
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SUBROUT livaS.. 


cee ce ee ee ce ce ce ee ee ee es ee ee ee ee es ee ee ee ce ee ee ee es es ee ee ee es ee ee ee ee ee es ee ee es ee ee es ee ee ee ee eee ee ee eee ee ee ee es 
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<r ee ee ee ee es ee ee ees ee es ee es ee es ee ee ee es ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee es ee ee ee ee ee es ee ee ee ee ee eee eee eee ee 
eee ee ee ee es ee ee ee ee es ee ees ee ee ee ee ee es ee ee es ee ee ee ee ee ee ee es ee es ee ee ee ee ee es es ee es ee ee ee ee eee ieee ee ee ee ome oe oe oe eee oe eee ss ee 


REAL FUNCTION PNEW(P,Q,MM,NN,A,B,C) 


ers SUNGTDON INTEGRATES THES HQUATIONeWITH RESPECT TO TIME 
eo DP/DT= (-D(UP)/DX -D(VP)/DY + A*DQ/DK + B*DEL(P) +C*V ) $ss 
UTILIZING AN UPWIND DIFFERENCING SCHEME 


DIMENSION P(MM,NN) ,Q(MM,NN) 
COMMON /PN/ DT,1I,J,UF,UFA,UB,UBA,VF,VFA,VB,VBA,VPT,FV 
COMMON /PARM/ H,M,N,P 


Pl= (UF-UFA)*P(I+1,J) +(UF+UFA-UB+UBA)*P(I,J) -(UB+UBA)*P(I-1,J) 
P2= (VF-VFA)*P(I,J-1) +(VF+VFA-VB+VBA)*P(I,J) -(VB+VBA)*P(I,J+1) 
P3= O(I+1,J3) -Q(I-1,J) 


P4= P(I+1,J3) +P(I-1,J)+P(1,J-1)+P(1,J+1)-4.0E0*P(I,J) 

P5= 2.0E0*H*VPT - ail 
TEMP1=(DT/(2.0E0*H))*((-P1-P2) +A*P3+(2.0E0/H) *B*P4+C*P5) 
PNEW=P(I,J)+TEMP1 

_ RETURN . 

END 


a ee ce ee ee ee ee ce ce ce ce ce ce ce ce ce ee ce ce ce ee ce ce ce ec ce ce ce ce ce ce ce ce ce ce cs ce ce ce ce ee es ee ee ce ee ee ee ee ee ee ee ee ee 
ce ce ce ce ce ce ee et ee ee ee ee ee es ee es ee ee ee es ee ee es ee ee es es ee ee es es es cs ee ce ec ey es ce ce ce ce es ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee oe 


REAL FUNCTION UBS(XU,YU,MM,NN,X,Y,ZETA) 


THISwSHERCUTENE FENDS BOUNDRY VALUES OF VELOCITY JU 
UTILIZING THE BIOT-SAVART LAW 


(NOTE: TERMS IN EQUATION ARE EVALUATED CONSECUTIVELY IN EACH 


4] 


QUADRANT FROM ONE TO FOUR ) 
DIMENSION X(MM),Y(NN) ,ZETA(MM,NN) 
COMMON /PARM/ HM ilare 


UMIN=-1, 02-05 
Os)... 00 
J2=.0.0E0 
USs= 102020 
U4= 0.0E0 


SIGNP=1.0E0 
SLGNM=— 12 Os0 
DA= H*H 


3 == Sas COMPUTE U COMPONENT OF VELOCITY US eG eee ool oo 


BO 20 t= 
DO 20 J= 1,N 
aa QUADRANT ONE -- | 
RSQ1=((XU-K(I))**2 + (YU+Y(J))**2) 
Ul= Ul+ ((-¥(J)-¥U)/(2.0E0*PI*RSO1 ) ) *STGNM*ZETA(I,J)*DA 
-- QUADRANT TWO -- 
RSOZ= ((XUfX(1I))**2 +(YU+Y (J) )**2) 
U2= U2+ ((-¥(J)-YU)/(2.0E0*PI*RSQ2) )*SIGNP*ZETA(I,J)*DA 
-- QUADRANT THREE -- . 
RSO3= ((KU+X(I))**2 +(YU-¥(J))**2) 
U3= U3+((Y(J5)-YU)/(2.0EO*PI*RSQ3) ) *SIGNM*ZETA(I,J)*DA 
-- QUADRANT FOUR -- 
'RSO4= ((KU-X(I))**2 + (YU-¥(J))**2) 
20 U4 = U4 + ((¥(3)-YU) / (2. O0EOXPI*RSO4) ) *SIGNP*ZETA(I,J)*DA 
UTEMP= U1+U2 +U3+ U4 
IF (ABS(UTEMP).LT. UMIN) UTEMP=0.0E0 
UBS= UTEMP 
RETURN 
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AAA AA A 


END 


REAL FUNCTION VBS(XV,YV,MM,NN,X,Y,ZETA) 


Poms Je. OUliN eer bibs BOUNDRY VALUES OF VELOCITY V 
UDEEECINGSERE BIOT-SAVART LAW 


(NOTE : TERMS IN THE EQUATION ARE EVALUATED CONSECUTIVELY IN EACH 
QUADRANT FROM ONE TO FOUR) 


DIMENSION X(MM) ,Y(NN) ,ZETA(MM,NN) 
COMMON/PARM/ H,M,N,PI 


VMIN= 1.0E-08 


Vil= 0.050 
V2= 0.0E0 
V3= 0.0E0 
_ V4= 0.0E0 
SIGNP=1.0EO 
SIGNM=-1.0E0 
= DA= H*H 


—s--" Saves, ) ONLONENT OF VELOCITY US@ENG THE B-S LAW---- 


DO 20 I= 1,M 
DO 20 J= 1,N 
-- QUADRANT ONE -- 
RSQL=((KV-X(I))*#*2 + (¥V+Y¥(J))**2) 
Vl= V1+ ((XV-K(1I))/(2.EO*PI*RSOL) )*SIGNMAZETA(I,J)*DA 
-- QUADRANT TWO -- 
RSO2Z= ((KV+X(1I))**2 +(YV4tY (J) )**2) 
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V2= V2+ ((XV+K(I))/(2.O0EO*PI*RSQ2) )*SIGNP*ZETA(I,J)*DA 
C, -- QUADRANT THREE -- . 
RSQ3= ((XV+K(1))**2 +(YV-¥ (J) )**2) 
V3= V3+((XV+K(1))/(2.0E0*PI*RSQ3) )*SIGNM*ZETA(I,J)*DA 
C  -- QUADRANT FOUR -- 
RSQ4= ((KV-X(1))**2 + (YV-Y(J))**2) 
20 V4 = V4 + ((XV-X(I))/(2.0Z0*PI*RSO4) )*SIGNP*ZETA(I,J)*DA 
VIEMP = V1+V2 +v3+ v4 
IF(ABS(VTEMP).LT.VMIN) VTEMP=0.0£0 
VBS= VTEMP 
RETURN 
END 


SUBROUTINE COUT(X,Y,MM,NN) 


C eeeceeuannan oun wena aunnaeananaans ee oe we we we es ee we ww oe we me we ee ee ee ee ee 
REAL X(MM),Y(NN) ,H,PI 
COMMON/PARM/H,M,N,PI 
WRITE(6,50) 

WRITE(6,*) | ARRAY X! 
WRITE (6,*) : 
WRITE(6,*) (X(I), I=1,M) 
WRITE (6,*) 
WRITE(6,*) 
WRITE(6,*) | ARRAY Y 1 
WRITE (6,%*) 
WRITE(6,*) (Y(J), J=1,N) 
WRITE(6,*) 
WRITE(6,*) 

50  FORMAT(///) 
RETURN 


a4 


CSBSSSSeSSSSe SaaS SS SSS SSS SSS SS SS SSS SSS SSS SSS SS SS SS SS SSS SSS SSS SSS SSS SSS 
SUBROUTINE OUTPT(A,LBL,MM,NN, KAL) 

C 

C THIS SUBROUTINE PRINTS OUT AN (NXM) ARRAY WITH A LABEL 

Q IN A RECTANGULAR ( 2,M-2)X(2,N-2) GRID 

C 


DIMENSION 2(MM,NN) 
COMMON/PARM/ H,M,N,PI 
CHARACTER*8 LBL 
NM2= N-2 
MM2= M-2 
WRITE (6,50) 
WRITE (6,*) 
WRITE(6,*) ! OUTPUT OF ARRAY: ',LBL,' ITER NUM:',KAL 
WRITE(6,*) | 
DO 10 J= 2,NM2 
WRITE (6,*) 
WRITE(6,*) (A(I,J), I=2,MM2) 
10 CONTINUE | | 
50 FORMAT('1') 


RETURN 
END 
C 
pe = Sr = - 3282+ “Se - - 9 -- ne ---- ------ 
C 
SUBROUTINE CPLOT(RHO, RHOI,MM,NN,NSTEP,TM,MW,NW) 
COMMON WK(15000) 
COMMON /PARM/ H,M,N,PI 
REAL W(81,81) 
REAL RHO(MM,NN),RHOI(MM,NN) ,H,PI,BHOT(50,100) 
DATA ISCALE/4HSCAL/ 
C 


eee eeeeey OF DENSEZTY CHANGES AT EACH NODE (DRHO) ....... 
e 


MM@e M=Z 
NM2= N-2 
DO 10 I= 2,MM2 
DO 10 J= 2,NH2 
JJ= (NM24+1)-J 
ITI= I+(M-5) 
RHOT(I,J)=RHO(I,J)+RHOI(I,J) 
W(IL, JJ)= ieee 0s) 
10 CONTINUE 
DORZO i= 370iiz 
II= (M-1)-I 
DO 20 J= 2,NM2 
JJ= (NM2+1) -J 
Zo W(II,JJ)= RHOT(I,J) 


CALL OUTPT(RHOT,"RAICUE a) a Nea) 
acta CALL SUBROUTINE CONTR’ TO PLOT THESE PERTUEe TONG)... 


CALL RESET (3HALL) 
CALL PAGE(8.7,11.2) 
CALL PHYSOR(2.0,3.00) 
CALL AREA2D (5.0,6.0) 
CALL HEIGHT(.2 ) 
CALL INTAXS 
| CALL ‘XNAME ('X/BO!,4) 
CALL YNAME ('Y/BO!,4) 
CALL KTICKS(2) 
CALL YTICKS(2) 
CHUMNGRAF (-5.0,1.,SsemGm a5 08) 
CALL MIXALF(! INSTRU") 
CALL MESSAG('T(EH.5)*(EXHX) =$',100,1.5,6.5) 
CALL REALNO(TM,106, ‘ABUT! , 'ABUT') 
CALL FRAME 
CALL BCOMON (15000) 
CALL CONMIN(2.0) 
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CALL CONANG (20.) 

CALL CONMAK (W,MW,NW,ISCALE) 

CALL HEIGHT(.08) 

GAUMACONLIN (O,S360LID ,'HAEELS' ,2,5) 
CALL CONLIN (1,4HDASH, 8HNOLABELS ,1,3) 
CALL RASPLN(0.25) 

CALL CONTUR (2,6HLABELS, 4HDRAW) 

ALL HEIGHT(.2) 

CALL COMPIX 

CALL RESET('HEIGHT') 

CALL RESET('COMPLX’ ) 

CALL ENDPL (0) 


RETURN 
END’ 
C 
EM ee eo 8 eee eae eee eset sl. -soee-~------- 
C 
- SUBROUTINE ZPLOT(ZETA,MM,NN,TM,MW,NW) 
REAL ZETA(MM,NN) ,H, PI | 
COMMON WK(15000) 
COMMON /PARM/ H,M,N,PI 
DIMENSION W(81,81) 
DATA ISCALE/4HSCAL/ 
c 


wea. CREATE ARRAY FOR PLOTTING OF ‘ZETA IN ACTUAL COMPUTATIONAL AREA.. 
; | | 
| MMZg=9 M=Z 
NMZ=N-Z ‘ 
BOT Oe= 2 MMzZ 
II= I+(M-5) 
DO 1@avi= 2aNliz 
JJ= (NM2+1)-J 
10 Cie) =ZETACT, J) 
DO 20 I= 3,MMZ2 
II= (M-1)-I 


~ 
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eee? ¢ 


DO 20 J= 2,NM2 
JJ= (NM2+1) - J 
Will ,JJ) = -Ze nace 


RESET (3HALL) 
PAGE(8.7,11.2) _ 
PHYSOR(2.0,3.00) 


AREA2D (5.0,6.0) 

HEIGHT(.2 ) 

INTAXS 

XNAME ('X/BO',4) 

YNAME ('Y¥/BO',4) 

eT) 

SATE) S172), 

GRAF (-5.0,1: 5 loa Gee 


MIXALF (' INSTRU! ) 
MESSAG('T(EH.5)*(EXHK) =$',100,1.5,6.5) 
REALNO(TM,106,'ABUT', ‘ABUT! ) 
FRAME 

BCOMON (15000) 

CONANG (60.) 

CONMIN(2.0) | 

CONMAK (W,MW,NW, ISCALE) 

HEIGHT (.08) 

CONLIN (0,5HSOLID, 'LABELS' ,2,5) 
CONLIN (1,4HDASH, 'LABELS',1,3) - 
RASPLN(0.25) 

CONTUR (2,6HLABELS, 4HDRAW) 
HEIGHT ( .2) 

COMPLX 

RESET ('HEIGHT') 

RESET (! COMPLX'! ) 
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CALL ENDPL (0) 


RETURN 
END 
C 
é 
nae een en ee 2 ------------------ 
a 
SUBROUTINE DRPLOT(RHO,MM,NN,NSTEP, TM, MW,NW) 
COMMON WK (15000) 
- COMMON /PARM/ H,M,N,PI 
REAL W(81,81) 
REAL RHO(MM,NN),H,PI 
DATA ISCALE/4HSCAL/ 
€ 


C...CREATE ARRAY OF DENSITY CHANGES AT EACH NODE (DRHO) ....... 
€ 
MM2= M-2 
NM2= N-2 
DO 10 I= 2,MM2 
DO 10 J= 2,NM2 
JJ= (NM2+1)-J 
I= I+(M-5) 
W(II,JJ)= RHO(I,J) 
cae CONTINUE 
DO 20 I=.3,MM2 
- TI= (M-1)-1 
DO 20 J= 2,NM2 
JJ= (NM2+1) -J 


20 Meine =eRHOlT , a 
é 
@e..... CALL SUBROUTINE CONTR TO PLOT THESE PERTUBATIONS........... 
C 


CALL RESET (3HALL) 
GAmeMPAGE (827) 11.2) 
CALL PHYSOR(2.0,3.00) 
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CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


AREA2D (5.0,6.0) 

HEIGHT(.2 ) 

INTAXS 

XNAME ('X/BO!,4) 

YNAME ('Y/BO' ,4) 

RE LEKSIC2) 

YTICKS (2) 

GRAF (-5 071) 5. oan 
MIXALF (' INSTRU’ ) 
MESSAG('T(EH.5)*(EXHX) =$',100,1.5,6.5) 
REALNO(TM,106,'ABUT', 'ABUT' ) 
FRAME 

BCOMON (15000) . 

CONMIN (2.0) 

CONANG (20.) 

CONMAK (W,MW,NW, ISCALE) 

HEIGHT (.08) : 

CONLIN (0,5HSOLID,'LABELS' ,2,5) 
CONLIN (1,4HDASH, 8HNOLABELS, 1,3) 
RASPLN(0.25) 

CONTUR (2,6HLABELS , 4HDRAW) 
HEIGHT (.2) 

COMPLX 

RESET('HEIGHT') 

RESET (!COMPLX'! ) 

ENDPL (0) 


RETURN 


END 


oo) 


= coe oe ce ee ec ge es ce ce ee ce es ce ee ee ee ee es ec ce cs es cs cs cs cs cs cs cs cs es es es es es es se es es cs es es es ees ees es es ee ee ee 
oa oe ee ee ee ee ee ee eee ee es ee ee ee ee ee ee ee ee es ee es ee ee es ee es ee ee ee es ee ee ees es ee es es es es es es es es es es ee es ees es es es es es es ees es ee ee es ee ee ie 


SUBROUTINE CON(Z2,XRANGE, YRANGE ,NX,H,MM,NN,ZMAX,ZMIN) 


THIS SUBROUTINE GENERATES A CONTOUR PLOT OF ZETA IN A RECTANGUEE 
DOMAIN OF SIZE: YRANGE * XRANGE ON THE PRINTER. 


50 


CHARACTER*1 SYMBOL(8) ,GRAPH(i00) 
REAL LX,LY 

DIMENSION Z2(HM,NN) , VALUE (7) 
DamemertsOL / Kee Se tt et ty 
VALUE(1)= ZMIN 

VALUE(2)= ZMIN*.75 

VALUE(3)= ZMIN*.5 

VALUE(4)= 0.0 

VALUE(6)= ZMAX*.75 

VALUE(7)= ZMAX 

VALUE (5)=ZMAX* .45 

NY= 0.8*NX*YRANGE/XRANGE 

DELK= XRANGE/NX 

DELY= YRANGE/NY 

WRITE(6,50) 

WRITE(6,*) 'NX=',NX, 'NY=!,NY,'KRANGE=',XRANGE, 


x 'YRANGE=! , YRANGE 
WRITE (6, *) | 
DO 11 I= 1,NX 
il GRAPH(I)= 'T! 


WRITE(6,6) (GRAPH(I), I= 1,NX) 
DO 5 JSYMBL = 1,NY 
(JSYMBL-0.5)*DELY 
2 14Y/H 
LY= Y-(J-1)*H 
HMLY= H-LY 
DO 4 ISYMBL = 1,NX 
X= (ISYMBL-0.5)*DELX 
IT=1+X/H 
LX= X-(I-1)*H 
HMLX= H- LX 
Al= HMLX*HMLY 
A2= HMLK*LY 
A3= LA*LY 


Cy 
I I) 


(rs a oe @ eee 
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60 


A4= LX*HMLY 
Z2c= (AY¥Z2(1,J) +A2*22(1, Jee eeeneeeerer oe 
* +A4%Z2(I+1,J))/H**2 . 


DETERMINE THE VALUE OF NRANGE BASED ON 220.0). eee 
PRINT OUT CHARACTERS IN ARRAY CREPH BY (OW 


00 21,7 
IF (Z2C-VALUE(K)) 1, 1, 2 
NRANGE = K 
ae ato 2 
CONTINUE 
NRANGE= 8 
GRAPH( ISYMBL)= SYMBOL(NRANGE) 
CONTINUE 
WRITE(6,6) (GRAPH(I), I=1,NX) 
FORMAT(' ',10X,100A1 ) 
FORMAT ('1') 
DO 60 I= 1,NX 
GRAPH(I)= 'B! | 
.WRITE(6,6) (GRAPH({I) ,I=1,NX) . 
WRITE(6,*) 
WRITE(6,*) 'MAXIMUM VALUE OF NEGATINE VORTICITY = SYMBOL (X) = ' 
x | 2MIN | . 
WRITE(6,*) ' 75% OF MAXIMUN NEGATIVE VORTICITY 
WRITE(6,*) ' 50% OF MAXIMUM NEGATIVE VORTICITY 
WRITE(6,*) 
WRITE(6,*) ' MAXIMUM VALUE OF POSITIVE VORTICITY = SYMBOL (1) =! 
* | ZMAX 
WRITE(6,*) ' 75% OF MAXIMUM POSITIVE VORTICITY = SYMBOL (+)! 
RETURN | 
END 


SYMBOL (A) ' 
SYMBOL (-) ! 


OM 


SUBROUTINE ZPL(ZETA,MM,NN,XWIDTH, YHT, ITER) 
COMMON /PARM/ H,M,N,PI 
DIMENSION Z2(50,100),ZETA(MM,NN) 


mere non leeeees SEOReeLOITING OF 2ZelA IN ACTUAL COMPUTATIONAL AREA.. 


MiSs has 
NM3=N-3 
Z2MIN= 0.0 
ZMAX= 0.0 
DO 10 T= litt 
DO 10 J= 1,NM3 
Z2\ oe —= ZEW I+1,J+1) 
fee oN THE MAXIMUM POSITIVE VORTLCITY iN ARRAY “EIS 


TF(Z2(1,J) .GT. ZMAK) ZMAK=Z2(1,J) 
Bee Ihe MAXIMUM NEGATIVE VORTICITY IN ARRAY Z2(MM,NN)....... 


C 
IF( Z2(I,3) .LT. ZMIN) ZMIN= Z2(I,J) 
10 ° CONTINUE : 
IF (ZMAX .LT. 1.0E-6 ) ZMAX=1.0 
. CALL SUBROUTINE CON TO PLOT VORTICITY.....-.cccseecseeceucs 
C a " ee . oak — a : ‘ 
CALL CON(Z2,XWIDTH, YHT,100,H,MM,NN, ZMAX, ZMIN) 
RETURN 
END 
Cee eee em ene ew em ew ew meme eee we we ew ee ee ee ee ee eee ee ew ee ee ee ee eee ee ee ee ee eee ee eee 


SUBROUTINE VPLOT(U,V,X,¥,MM,NN,TM) 

COMMON /PARM/ H,M,N,PI 

COMPLEX ZV(50,100), CV(50,100) 

DIMENSION U(MM,NN), V(MM,NN), X(MM), Y(NN) 
MM2= M-2 

NM2= N-2 


Oa Cie @ (Se @ aaa ee 


MM3 =M - 3 
NM3 =N - 3 
C 
one FORM COMPLEX U, V, AND Z ARRAYS FOR QUADRANTS 3 AND 4 ...... 
é 
DO 10 I= 2,MM2 
DO 10 J= 2,NM2 
Tl ager 
JJ= (NM2+1) -J 
ZV(II,JJ)= CMPLX(X(I),¥(J)) 
10 CV(II,JJ)= CMPLX(U(I,J),V(I,J)) 
C . 
C | 
Coane PLOT VECTORS USING SUBROUTINE V2PLOT ..........ecceeees 
CALL V2PLOT(ZV,CV,MM3,NM3,7M) 
C 
RETURN 
END 
C 
Ca owen cnececsccceseeleeccce see seeeee se ee 
C i 
SUBROUTINE V2PLOT(Z,C,M,N,TM) 
RAEKKAARAKKARKAKRAKRKKAKRRRAARAAAKAKAAKAKAKAAKKKAAAKAKNARAKRKRRKRRARARKRR 
* SUBROUTINE V2PLOT TO PLOT MKN COMPLX VELOCIT FIELD C(M,N) * 
* OVER A COMPLX DOMAIN Z(M,N) _ 
x x 
KRAKREKKKRKRKARRKRKKKKKKAKAKAKARAARKRAKKRARKKAKKRRKRRRR RAR RRR RRA KRARKRKRKKARKARKX 
COMPLEX Z(50,100),C(50,100) , CMAX 
REAL X(100),¥(100) ,¥X(100) ,¥N(100) ,XX(100), 
* ¥N(100) ,X2(100),¥2(100) 
INTEGER PMAX, PMIN 
EXTERNAL AMAX 
EXTERNAL AMIN 
C 


DO 20 J =1,N 
DO 10 I =1,M_ 
X(I) =ABS(REAL(C(I,J))) 
¥(I) =ABS(AIMAG(C(I,J))) 
CONTINUE 
CALL AMAX(Y,M,YMAX, PMAX) 
CALL AMAX(X,M, XMAX, PMAX) 
YX(J) = YMAX 
XX(J) = XMAX 
CALL AMIN(Y,M,YMIN, PMIN) 
CALL AMIN(X,M,XMIN, PMIN) 


YN(J) =YMIN 
XN(J) =XMIN 
CONTINUE 


- CALL AMAX(YX,N, YMAX, PMAX) 

CALL AMIN(YN,N,YMIN, PMIN) 
CALL AMAX(XX,N,XMAX, PMAX) 

CALL AMIN(XN,N,XMIN, PMIN) 


CMAX = CMPLX(XMAX, YMAX) 
UMAX =SQRT((XMAX)**2+(YMAX)**2) 
DO 24 I =1,M | 
DO 22 J=i,N 
S@o ys) neers) /UAK 
_ CONTINUE 
CONTINUE 


DO 40 J =1,N 
Bees0 I =1.,M 
Mee =RbAL(ZCL 0) >) 
V(l) =ADMEG(Z(1,3)) 
CONTINUE 
CALL AMAX(Y,M, YMAX, PMAX) 
CALL AMAX(X,M, XMAX, PMAX) 
YX(J) = YMAX 


Car 
Car 


40 


210 


XX(J) = KMAX 
CALL AMIN(Y,M, YMIN, PMIN) 
CALL AMIN(X,M,XMIN, PMIN) 
YN(J) =YMIN 
XN(J) =XMIN 
CONTINUE 
CALL AMAX(YX,N,Y¥MAX,PMAX) 
CALL AMIN(YN,N, YMIN,PHIN) 
CALL AMAX(XX,N,XMAX, PMAX) 
CALL AMIN(XN,N,XMIN, PMIN) 


CALL RESET(3HALL) 
CALL PHGHCS.7 11.2) 
CALL PHYSOR(2.0,3.0) 
CALL AREA2D(5.,6.) 
CALL HEIGHT(.2) 

CALL INTAXS 

CALL XNAME('X/BO',4) 
CALL YNAME('Y/BO' ,4) 
GALL XURGRS (2) 

CALL YTICKS(2) 

CALL GRAE(0.0,0.5 Sec) leon ino momen 


CALL MIXALF(! INSTRU! ) 


CALL MESSAG('T(EH.5)*(EXHA) wae 10C piece 
CALL REALNO(TM,106,'ABUT', 'ABUT' ) 


DO 60 J=1,N 
DO 50 I=1,M 
X(I) = REAL(Z(I,J)) 
X2(1) = X(I)+REAL(C(I,J)) 
Y(I)= AIMAG(Z(I,J)) 
Y2(1) = Y(I)+AIMAG(C(I,J)) 
CONTINUE 
CALL MARKER (3) 
Chil sGerie(. eur 


56 


CALL CURVE(X,Y,M,-1) 
DO 54 I =1,M 
a= aa) 
Ma 2) 
Y11=Y(I) 
Y12=Y2(1) 
IF(X12.£0.0..AND.Y12.£0.0) GO TO 54 
enews Sema) Yi) )X12,Y12,1121) 
54 CONTINUE 
60 CONTINUE 
CALL RESET(!HEIGHT') 
CALL ENDPL(0) 


RETURN 
END 
C 
C 
SUBROUTINE AMAX(Y,N,YMAX,PMAX) . 
: | 
G ee ERA RARER AA A RRA RRA RRERERARRRRR 
C * SUBROUTINE TO COMPUTE THE LARGE ELELEMENT 
C * IN A GIVEN ROW OR COLUMN IN ARRAY 'A!. 
e ee RRR AAA REAR RRKKKK RRR A 
c 
@ AX VARIABLE DECLARATION *** 
C REAL AMAX 
' INTEGER PMAX 

DIMENSION Y(N) 
E 

YMAX=Y (1) 

DemsOOmt—72 hi 

IF(Y(I).LE.YMAX) GO TO 5000 
YMAK=Y (TI) 
PMAX=I 

5000 _ CONTINUE 
5100 CONTINUE 


e MAKXROW=POS 


Na AAA A AN A 


RETURN 
END 
C 
C 
C 
SUBROUTINE AMIN(Y,N, YMIN, PMIN) 
KRM KKK NK KR ARKMRRKA MK KAKA KAMKKMAKRAAKAARE 
* SUBROUTINE TO COMPUTE THE SMALL ELEMENT 
x IN A GIVEN ROW OR COLUMN IN ARRAY 'A'. 
FoF KIKI KIRK RRR KR RERRKKRRRKRRRRRRKRREKKKKRRARRARRR 
*k*k VARIABLE DECLARATION *** 
REAL AMIN 
INTEGER PMIN 
DIMENSION Y(N) 
; ; 
YMIN=Y(1) 
DO 6100 I=2,N 
LF(Y¥(I).GE.YMIN) GO TO 6000 
YMIN=Y¥ (1) 
PMIN=1 
6000 CONTINUE 
6100 CONTINUE 
RETURN 
END 
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Figure B.9b Vorticity Field: SP = 0.00. 
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Figure B.13b  Vorticity Fieid: SP = 0.00. 
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Figure B.14 Experimental and Numerical Results for Various Initial Core Radii. 
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Figure B.15a Constant Density Contours: SP = 1.00. 
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Figure B.15d Vorticity Field: SP = 1.00. 
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Figure B.l6a Constant Density Contours: SP = 1.00. 


20 


T" - 


UO. 54580 


AAAAAAAAAAAAAAAAAAAAARAAAAAAAAARAAAA AAA AA AAA AAA AAI AVIV OTITIS 
A AAANAAAAAAAAAAAAAARAAAAARAAAAARAANAAA AAA AD AAA AA AAAI IATA hh 
AAKAAARARAAAAAAAAAAARAARRAAAARARA AA AANAA AAA AA NAAAAAADADI TTI 
AAAAAAAAAAAARARARAAAAARAAANAAAAAAANANA AAA NAA AA AAA AAIDAAA TOT IT 
AA AAAAAAAAARAAAARAAAAAAAAAAAAAAARAAAANA ANA AA AAA AAA AA AIVAA TTI Nan 
NA AAAAAAAAAAAAAAAAAAAAARARRAARAARAAAANAN AAA AAA AI AA AAAI PTO aa 
AAA AAAAAKAAAAAARAAAARAAAAAAARAARARARARAAAA AANA AA AAAIIAD IIT IIs 
AAAAAAAAAAARAAAAARAARAAAARAAAAAARANANAAAANA AY AAA AAA ANADAIDI TT WD 
AAAAAAAAAAAAAAAAAARRARARAAAAAAAAAAAAANAR AAA AAAA IAA AA AIDA IIT IN WE 
AAAAAAAAAAAAAARAAAAAAAAAAAAAAAAAAAARAAAN AA ADAYA AAA AAADATI IAAT DN 
AAAAAAAAAAAAAAAAARAAAAAAAAAAAARAAAARAAAAAN AA AAA AAA AA AAD IAAI ANAS 
AAAAAAAAAAAAAAAAAAARAAAAARARARAAARAARAAAAANAA AA AAA AAA ADDIS P FIA 
ANAAAAAAAAAAAAAAAARARARAAAAAARAAARAAAAARA AAAS AAJA AAA AAD ID ITIP AANA 
AAAAAAAAAAAAARAAAAAARAAAAAAAARAAARARAARRARNAAAAA AS AA AA FDA AIA 2 AA WNT 
AAAAARAAAARAAAAAARRAAARAAAANARANAAAAAAARAANAAAAAA AAA FADIA TTF WA 
AAAAAAAANAAAAARARAAAAAAAAAAARAAARAAAARAANANAANAA DAD AI AAAD IATA 72 OWI 
AAAAAAAANARAAAAARANRANAAAAAARAAAAAAARAAAAAAAAAAAAA AAA AA AFI TIE OWI 
ANAAAAARAARAANAARAAAAAAARRAARAANRAAAAAAARAAAAAARAAAAAAAAAAD ITA OATS 
NAAKAAAAAAAAAAAARRRRARRAAARARARAAARRAAAAAARAARAAAA AAD AA AAAI II ISO ANNA 
AAAAAAAAAAAAARAARARARARARARARAARARARAAARAAAAAAANAAAA AAA AAA TI 7772 > OANA 
NAAAAAARAAANAAAAAAAAAAAANRARAAARARAAAAARAARAAARAAAAAAA AAAI 99772 OA VY 
AAAAAAAAAAAAAAARARAAAAAAAAAAAANARARAAAARAAAAAAAAAAAAS4A A AAAI FIFO ODOANMNAANTN YY VV VV 
AAAAARAAAANARAAAAAAAAAARAARAARARAAARARAARARAARAAAAANAAA 444 AFD SANNA 
AAAANAAARAANAAARRARAAAAAARARRARRAAAARRAARAANARAAAAAAAA AAA AA 44-47 77 ANN AV 
AAAANAARAAAAARARRRNAAARRRAAAAARARAAARARRAARARAARAAAAAAAA4 49-97 7 AAA 
AAAAAARAAARARARRARRRABIARRRRA AAARBRARAR AAA BAARANNARAAAAAAAAA4 44773 Aan YY VV VY 
AAAAARRAAARAAARARRARARR DRI MARRIARBAI ARR AARRAARANAAAAAAAAAAAA 44 44-775 >3 nas aqayyyVVYVVVV 
AA AARAAAARRNARRRRRRIRDN AAD DRADER ARARIARADDAARAAAAAA AMAA A A Areas g guy ywYVVVVVY 
DAKAR AAR AARKRDEDARERE EDD DRED DDD DARA RDRARRA DAA RAA AAA ADD Aeros yyy VYWEV YY 
AN AAR AA AAARKSDNDREDDDEREEEEDDR DEDEDE DERRIDA ABABA BAR MARAE ET terra yy VRE 
AAA AAR ARK RARARRESREREDDRERER ERED DDR ERD EDR DARARD ARR AAA AM PRAA TES VV EY 
AAA ARIANA SREREDEEREEREDEDEDEDENDREDEDRDDI DD BA AA AAAS WA ZN \ a aaa 
AAAAARRREEDEAEDERELEDELEDEEEDEEDEEEDE REDD EDDA A AAA AE F ow Winaaaaad 
AARARDRDDDDENETTEEDTTTLTEDEPEETOEPLLTDEDEDPDD ND PANNE A A IN f ay Vi pppe bere 
f- \\ 






AA ARRREREREE REPRE EEE PE PEE ETT TET HEE EPEHEBEPPWHH HH PH 





ANKARS PREPEPEE TPT TT PEE TELE TET EE EE EEEE EEE EP EREERE RD \ " a pVeREELEL 
AAARBESYEYD EEE EERE EEE IEEE EEE EEE TEEPE PEEP P PRIN : LW, JU bbbhhhse 
AD ARRRELE STRAT NTT STE OE EE ETT UE EEE PEE Se eg & jt bebe be bbe bok 
AREEYTTT Be WWE ene he 

\ RUNV SSS SS SSO cdc cdecaccccccce Sas is 5 OTTO ae: 


“~ 


» 


U a a ee Ol ee So ee lee 


265 
K/ BO 


: 


EE ee 


Figure B.16b Velocity Field: SP = 1.00. 
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Figure B.16d Vorticity Field: SP = 1.00. 
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Figure B.l7a Constant Density Contours: SP = 1.00. 
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Figure B.l7c Density Perturbatton Contours: SP = 1.00. 
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Figure B.l8a Constant Density Contours: SP = 1.00. 
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Figure B.18d Vorticity Field: SP = 1.00. 
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Figure B.19a Constant Density Contours: SP = 1.00. 
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Figure B.21a Constant Density Contours: SP = 1.00. 
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Figure B.2Ic Density Perturbation Contours: SP = 1.00. 
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Figure B.22a Constant Density Contours: SP = 1.00. 
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Figure B.23a Constant Density Contours: SP = 1.00. 
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Figure B.24a Constant Density Contours: SP = 1.00. 
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Figure B.24b Velocity Field: SP = 1.00. 
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Density Perturbation Contours: SP = 1.00. 
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Figure B.24d_ Vorticity Field: SP = 1.00. 
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Figure B.25a ‘Constant Density Contours: SP = 1.00. 
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Figure B.25b Velocity Field: SP = 1.00. 
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Figure B.25¢ Density Perturbation Contours: SP = 1.00. 
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Figure B.25d Vorticity Field: SP = 1.00. 
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Figure B.26a Constant Density Contours: SP = 1.00. 
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Figure B.26b Velocity Field: SP = 1.00. 
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Figure B.26d Vorticity Field: SP = 1.00. 
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Figure B.27a Constant Density Contours: SP = 1.00. 
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Figure B.27b Velocity Field: SP = 1.00. 
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Density Perturbation Contours: SP = 1.00. 
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Figure B.27d Vorticity Field: SP = 1.00. 
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Figure B.28b Velocity Field: SP = 1.00. 
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Figure B.28d Vorticity Field: SP = 1.00. 
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Figure B.29a Constant Density Contourssor = 1°00. 
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Figure B.37a Constant Densitv Contours: SP = 0.50. 
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Figure B.38a Constant Density Contours: SP = 0.50. 
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Density Perturbation Contours: SP = 0.50. 
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Constant Density Contours: SP = 0.50. 
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Figure B.40b Velocity Field: SP = 0.50. 
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Figure B.40c Density Perturbation Contours: SP = 0.50. 
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Figure B.4lc Density Perturbation Contours: SP = 0.50. 
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Figure B.42a Constant Density Contours: SP = 0.50. 
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Figure B.42b Velocity Field: SP = 0.50. 
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Figure B.45b Velocity Field: SP = 0.50. 
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Figure B.45d_ Vorticity Field: SP = 0.50. 


Y/BO 


T” = 0.59520 





0.0016 ———__ 0.0016 


-_—_=_ = ee ee ws eee ewe ee eee ee ee See ee ee eS eee ee ee ee eee eee eee eee see i i ll 


0.002 eee 


I ee ee ee ee ee ee ee ee ee ee ee ee ee ee 





X/ BO 
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Figure B.46b Velocity Field: SP = 0.50. 
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Figure B.46c Density Perturbation Contours: SP = 0.50. 
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Figure B.46d Vorticity Field: SP = 0.50. 
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Figure B.47a Constant Density Contours: SP = 0.50. 
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Figure B.47d Vorticity Field: SP = 0.50. 
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Figure B.48c Density Perturbation Contours: SP = 0.50. 
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Figure B4Sd Vorticity Field: SP = 0.50. 


Vass 


Y/BO 


10 





—_ = Se = . 
—_ 7 
—_ 
- ~~ «as 
a Ad —_— 
== —~— = =—= 
—_ = — a= 
—_—__e_we es = —=—— ew ie i 


=~ wee wesw wes eS ’ = ie i i ee ee ee Ce 


0.0016 =. ee CU 


eee eee eee ee i i se i i o_o 


0.0020 = ae 


em mm ee eel lie 
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